set more off
program drop _all
program define REpbias, rclass
version 13.1

syntax, alpha(real) studies(integer) select(real) obs(integer)

// Remove existing variables
drop _all

//We first create the matrix to store the results of each study
matrix A = J(`studies',1,0)
matrix B = J(`studies',1,0)
matrix C = J(`studies',1,0)
		
forvalues i = 1/`studies' {
clear
// STEP ONE: Create the data for each study and estimate an effect
set obs `obs'
generate x = rnormal()
// Note that each "study" has the same number of observations (100)
// but differ in the variance of their respective error terms.  This
// causes the estimate of the effect to be estimated with varying degrees 
// of precision
scalar lambdai = 0.5+30*runiform()
generate e = lambdai*rnormal()
scalar alphai = `alpha' + rnormal()
generate y = 1 + alphai*x + e
quietly regress y x 
scalar coef = _b[x]
scalar secoef = _se[x]
scalar tcoef = coef/secoef
	
//if abs(tcoef) < 2 {
// An alternative criterion is that results that show a negative effect
// have a harder time getting published.  To study that case, substitute
// the line below for the line above
if coef < 0 {
scalar dummy = cond(runiform()<`select',1,.)
// The statement above creates a dummy variable that randomly selects which
// studies will get "published" if they fail to meet the "publication criterion"
// either (i) abs(tcoef) >= 2 or (ii) if coef >= 0.  Studies that are not "published"
// receive missing values and thus are not included in the "meta-analysis."
 scalar coef = dummy*coef
scalar secoef = dummy*secoef
scalar tcoef = dummy*tcoef
}



matrix A[`i',1] = coef
matrix B[`i',1] = secoef
matrix C[`i',1] = tcoef
}
		
// The next set of commands moves the data out of matrices and reformats them as 
// standard Stata data series.  We have now completed generating our individual
// studies and we now move into the "meta-analysis" stage.
matrix bob = A,B,C	
svmat bob
rename bob1 effect 
rename bob2 seeffect 
rename bob3 teffect 
generate pet = (1/seeffect)
	
// This estimate produces the PET version of the "publication bias"-
// corrected effect estimate
regress teffect pet, vce(robust)
return scalar effect_PET = _b[pet]
scalar effect_PET = _b[pet]
test _b[_cons] = 0
return scalar pvalue_FAT = r(p)
scalar pvalue_FAT = r(p)
test pet = `alpha'
return scalar pvalue_PET = r(p)
scalar pvalue_PET = r(p)
test pet = 0
return scalar pvalue_PETFPP = r(p)
scalar pvalue_PETFPP = r(p)

// This estimate produces the PEESE version of the "publication bias"-
// corrected effect estimate
regress teffect seeffect pet, noc vce(robust)
return scalar effect_PEESE = _b[pet]
scalar effect_PEESE = _b[pet]
test pet = `alpha'
return scalar pvalue_PEESE = r(p)
scalar pvalue_PEESE = r(p)
	
return scalar effect_FPP = effect_PET
return scalar pvalue_FPP = pvalue_PET

if pvalue_PETFPP < 0.05 {
return scalar effect_FPP = effect_PEESE
return scalar pvalue_FPP = pvalue_PEESE
}	

// This last command keeps track of how many studies are in our "meta-analysis"
return scalar N = e(N)

end
